Comparison of qPCR and metagenomic sequencing methods for quantifying antibiotic resistance genes in wastewater

Surveillance methods of circulating antibiotic resistance genes (ARGs) are of utmost importance in order to tackle what has been described as one of the greatest threats to humanity in the 21st century. In order to be effective, these methods have to be accurate, quickly deployable, and scalable. In this study, we compare metagenomic shotgun sequencing (TruSeq DNA sequencing) of wastewater samples with a state-of-the-art PCR-based method (Resistomap HT-qPCR) on four wastewater samples that were taken from hospital, industrial, urban and rural areas. ARGs that confer resistance to 11 antibiotic classes have been identified in these wastewater samples using both methods, with the most abundant observed classes of ARGs conferring resistance to aminoglycoside, multidrug-resistance (MDR), macrolide-lincosamide-streptogramin B (MLSB), tetracycline and beta-lactams. In comparing the methods, we observed a strong correlation of relative abundance of ARGs obtained by the two tested methods for the majority of antibiotic classes. Finally, we investigated the source of discrepancies in the results obtained by the two methods. This analysis revealed that false negatives were more likely to occur in qPCR due to mutated primer target sites, whereas ARGs with incomplete or low coverage were not detected by the sequencing method due to the parameters set in the bioinformatics pipeline. Indeed, despite the good correlation between the methods, each has its advantages and disadvantages which are also discussed here. By using both methods together, a more robust ARG surveillance program can be established. Overall, the work described here can aid wastewater treatment plants that plan on implementing an ARG surveillance program.


Introduction
Since their first use almost 100 years ago, antibiotics have been broadly employed in various fields such as agriculture, aquaculture, animal husbandry and medicine [1,2].Unfortunately, this has led to a global increase in antibiotic resistance genes (ARGs), collectively known as the "resistome" [3,4].Due to the misuse of antimicrobial drugs and compounds, ARGs and antibiotic resistant bacteria (ARB) have spread extensively in different environments [5].The World Health Organization (WHO) has classified ARB and ARGs as one of the top issues threatening the public health in the 21 st century [6].Alarmingly, it is predicted that by 2050, annual deaths as a result of infections with antibiotic resistant pathogens will surpass 10 million [7].In the most positive scenario where antimicrobial resistance is limited, it is projected that the global economy will experience losses of more than $1 trillion per year after 2030, and this number will increase to $2 trillion per year by 2050 [8].
ARGs could be circulated among bacterial cells either by horizontal gene transfer through mobile genetic elements (MGEs), such as plasmids, transposons and integrons [9,10] or vertical transmission via bacterial propagation [11,12].The selective spread of resistance genes occurs when bacterial cells are exposed to antibiotics [13].A lot of this exposure results from antibiotics that are discharged to the environment from animal and human wastes, industrial plants and antibiotics residuals in wastewater treatment plants [14][15][16].
Significant amounts of antibiotics have been discovered in water supplies, and numerous studies have shown the prevalence of ARGs in wastewater treatment plants (WWTPs) [17][18][19][20].Since the human body cannot metabolize around 60-85% of antibiotics, they are excreted into sewage wastewater [21].However, traditional biological treatment approaches in WWTPs are insufficient in fully removing ARGs and antibiotics as they are mainly designed for the elimination of solids and other organic matter [22][23][24].Importantly, since wastewater from different parts of a city all end up at large WWTPs, these WWTPs represent an opportunity to set up an effective ARG surveillance system.However, the success of this system is dependent on the ability to accurately characterize the resistome in these waters.
Recent advances in next-generation sequencing technologies and accompanying computational methods have become increasingly important for identifying ARGs.In read-based methods, ARGs can be detected by aligning whole-metagenome sequencing (WMS) reads to a reference database using assembly algorithms.A number of those tools have been made available in recent years [25].These technologies have provided additional benefits to the traditional PCR methods that are restricted by the number of primer pairs used, in comparison to the more extensive number of ARG genes available through resources such as the Comprehensive Antibiotic Resistance Database (CARD) [26], which we use in this study.In addition, sequencing-based methods enable downstream analyses that can help identify new mutations, phylogenetic diversity, genetic epidemiology, and potential loss/gain of function in ARGs can be inferred.
Sequencing-based methods have recently been recognized for their significant potential in studying the resistome [27].However, few efforts have been made to critically evaluate the advantages and disadvantages of PCR and sequencing methods for ARG identification in wastewater samples.To address this, we analyzed the ARGs in four wastewater samples from different locations in the Emirate of Abu Dhabi, UAE.We then compared the resulting composition of the ARGs detected by the Resistomap (WaferGen HT-qPCR) array to TruSeq DNA PCR-free sequences analyzed using the Resistance Gene Identifier (RGI) method offered by CARD.The RGI method uses protein data based on homology to predict resistomes.A qualitative and quantitative evaluation of the resulting genomic data was carried out, and a comparison between both methods' validity and sensitivity is discussed here.The comparison highlights the advantages of sequencing-based methods, which are able to cover the entire sequence of the gene by the assembly of reads and report mutations inside the genes.In contrast, the Resistomap qPCR array reports only the presence or absence of the limited genes they have designed primers for.
This study aims at providing recommendations on which genes to test for and what quality control filtering thresholds to use when combining both methods in order to achieve a balance between sensitivity and affordability.These recommendations will guide the implementation of ARG surveillance programs.

Sample collection
Four wastewater samples were collected from different location types in Abu Dhabi City, United Arab Emirates.Sample 1 (S1) was collected from an urban residential site, Sample 2 (S2) from an industrial site, Sample 3 (S3) from a healthcare facility, and Sample 4 (S4) was collected from a rural sewage treatment plant influent.No permits were required for sampling municipal wastewater in the Emirate of Abu Dhabi, as all the sampling was performed by employees of the Abu Dhabi Department of Energy and the Emirates Bio-research Center at the Ministry of Interior.
Composite wastewater samples were collected using an auto-sampler over a 24-hour period.250 ml of wastewater was transferred from the auto-sampler to sterile polypropylene bottles (ISOLAB, GmbH) on ice to the laboratory.Once received at the laboratory, samples were stored at 4˚C and gDNA was extracted within 24 hours.Municipal wastewater samples have been collected from all over the United Arab Emirates by our laboratory since April 2020 [28].Wastewater characteristics have been pretty consistent throughout this period with similar readings observed from wastewaters sampled from manholes or wastewater treatment plants.Municipal wastewater samples reported average values for NH 3 -N of 29.3 mg/L, PO 4 3-P of 14 mg/L, COD of 418 mg/L, conductivity of 3160 mS/cm and pH of 7.5 [29].

Wastewater gDNA extraction
Wastewater gDNA extraction was performed using the PowerSoil 1 DNA Isolation kit (MoBio Laboratories Inc., Carlsbad, CA, USA) following the instructions supplied with the kit with some minor modifications as described previously [29].Briefly, conical tubes, each containing 50 mL of wastewater sample, were mixed vigorously and then centrifuged (5,000 xg) for 15 min to obtain bacterial cell pellets.The pellet was suspended in the PowerBead buffer tubes (one of the DNA extraction kit components).At this point, the manufacturer protocol was followed.The concentration and purity of extracted gDNA was determined using a Nano-Drop™2000 spectrophotometer (Thermo Fisher Scientific™) and gDNA was stored at -20˚C until it was further processed.

SmartChip HT-qPCR assay
Wastewater gDNA samples were sent to Resistomap laboratory (Helsinki, Finland) to perform the HT-qPCR-based ARG assay.All runs were carried out using all 384 primer sets that are utilized by Resistomap [30].Of these 384 primer sets, 268 are used to detect ARGs that confer resistance to nine classes of antibiotics including aminoglycoside, beta-lactam, phenicol, macrolide-lincosamide-streptogramin B (MLSB), quinolone, sulfonamide, tetracycline, trimethoprim, and vancomycin.In addition, 38 primer sets are used to detect ARGs that confer multidrug-resistance (MDR), 17 primer sets to detect genes that confer resistance to antibacterial agents, and 52 primer sets for genes linked to mobile genetic elements (MGEs) such as transposons, plasmids integrons and insertion sequences (IS).Finally, to gain insight into the bacterial community size and structure, 1 primer set was used to detect total 16S rRNA genes and 8 primer sets for determining bacterial taxonomy (Firmicutes, Bacteroidetes, Acinetobacter baumannii, Pseudomonas aeruginosa, Campylobacter, Enterococci, Klebsiella pneumoniae and Staphylococci).A list of the primers used in this study was extracted from Resistomap's publication [30] and are provided here as supplementary data (S1 Table ).
The HT-qPCR reactions were carried out in the SmartChip Real-time PCR system (Wafer-Gen Biosystem, Freemont, CA, USA).This system uses 5184 well chips (100 nL each), which are auto-filled by the SmartChip Multisample Nanodispenser using amplification and detection parameters previously described [31].The results were analysed by SmartChip qPCR Software (v 2.7.0.1).The detection limit was set at a threshold cycle (CT) of 27 [32][33][34].The specificity of the primer sets was monitored through the melting curve of each PCR product.PCR efficiency values ranged between 1.60 and 2.35, with 95.5% of the values falling between 1.8 and 2.1.The mean, median and standard deviation of these values were 1.92, 1,94, and 0.12 respectively.The abundance of the discovered genes normalized to the total 16S rRNA genes per gDNA sample was calculated via relative quantification using the 2 -ΔCT method (Eq 1) [35].

TruSeq DNA sequencing
Wastewater gDNA samples were shipped to Macrogen Inc. (Seoul, South Korea) for library preparation and sequencing.The quantification of different samples was performed by the picogreen approach via Victor 3 fluorometry, and their quality was evaluated using gel electrophoresis.The PCR-free Illumina TruSeq DNA libraries were prepared according to the manufacturer's instructions (Illumina, San Diego, CA, USA).The libraries were adjusted to 4 nM and sequenced using the Illumina HiSeq sequencing system.

Shotgun sequencing and bioinformatics approach
FastQ files were analysed using the standard quality control tools FastQC and Qualimap [36].Subsequently, RGI-Resistance Gene Identifier (version 5.2.1), was used to map reads from all four samples against the Comprehensive CARD, version 3.2.2 using the command rgi bwt.An example of the code as used for sample S1 is provided below: rgi bwt-read_one S1_1.fastq.gz-read_twoS1_2.fastq.gz-alignerbowtie2-output_file AFY01-threads 30-local It invokes 'bowtie2'version 2.4.5 to perform reference-based assembly, using CARD as a reference which includes a total of 4891 known ARGs.The exact command issued by rgi is: bowtie2-align-s-wrapper basic-0-quiet-very-sensitive-local-x bowtie2 -S S1.temp.sam-1 S1_1.fastq.gz-2 S1_2.fastq.gzThis process performs conventional Burrow-Wheeler Transform mapping and alignment of reads, which yields an output in bam-file format.Subsequently, resulting tabular gene match reports are combined and visualized using Python/Jupyter (see GitHub repository rgivisual.ipynb).The gene map reports include for each gene the number of matched reads, the percentage to what extent the gene is covered by at least one read, the mechanism of resistance, potentially known hosts, the drug class, the resistance categorization and an ontology identifier for the Antibiotic Resistance Ontology (ARO).

Quality control and rarefaction
To determine whether a gene is detected or not, the percentage coverage of the gene length by reads is used.The rationale is to be able to filter by what fraction a gene has been detected.We also calculate fold-coverage, in order to get a sense of abundance as shown in Eq 2: We adapted the concept of rarefaction as commonly applied in microbial ecology, but instead of assessing the completeness of samples species, we quantified how comprehensive the sampling of ARGs and drug classes are.Rarefaction was conducted by randomly subsampling reads from the complete sample S1.This procedure yielded forward and reverse fastq files for subsets of 10%, 20%. . .90% of reads, with 3 replicas respectively, while ensuring that forward and reverse reads are subsampled as pairs (rarify.py in the specified GitGub repository).Subsequently, the shell script rgibwt2_rarify.shdetails the RGI runs for each individual dataset using parallelization with the same settings as the complete dataset.Lastly, the Jupyter Notebook rarefactionViz.ipynbvisualizes rarefaction curves by combining all RGI gene mapping reports.

Taxonomic classification
The Kraken2 tool [37] is a software program used in the metagenomic analysis for taxonomic classification of DNA sequences obtained from mixed microbial communities.Kraken2 is used for assigning taxonomic labels to the sequenced reads for the four samples using the following command: kraken2-threads 28-db $DB-output output.txt-reportreport.txt-paired$1.fastq.gz$2.fastq.gzwhere DB refers to the Standard-16 collection containing archaea, bacteria, viral, plasmid, human and UniVec Core capped at 16 GB, and $1 and $2 refer to forward and reverse reads for a given sample.

Post-processing of results for comparability
The comparison of reported drug classes between the PCR Amplification method and the CARD/RGI approach is hampered by the different naming conventions and drug class definitions used.In addition, drug classes are reported on very different levels of classification, however, due to CARD's ARO annotation, it is possible to infer high level classifications.We developed an Ontology-term based mapping algorithm, which establishes a link between 11 high-level Resistomap drug classes and the equivalent terms in the ARO annotation.For each Resistomap drug class, the algorithm then aggregates reported CARD genes that are subsumed under the corresponding ARO class (direct matches or direct ontology children).This was done with the help of the OBO parser from the python library GOATOOLS [38], operating on the aro.obo file as retrieved from OBO-Foundry (see GitHub repository linkResistome2Card-DrugClasses.py).For example: The Resistomap drug class "MLSB" is mapped to the three ARO classes Macrolide antibiotic (ARO:0000000), Lincosamide antibiotic (ARO:0000017) and Streptogramin (ARO:0000026).
Likewise, gene naming conventions in the two datasets differed, necessitating establishing equivalent genes reported.This was accomplished by matching Resistomap's PCR primer pairs against the CARD database (script scanPrimers.py).In addition, when primer pairs matched multiple targets, we analysed, how related these targets are with the help of the ARO graph: we extract the percentage of hits that have the same (grand-)parent node as well as the maximal distance between the hits.For the most different hit pair, we report the most recent common ancestor (using Python/networkx method all_pairs_lowest_common_ancestor, see summarizeLinksARO.py in GitHub repository).
Pearson Correlation was used to determine the degree of correlation between the AMR genes identified by both methods.Pearson correlation coefficient was calculated for both the unique primer and multi-target primer matches.Furthermore, the breakdown of the correlation by antibiotic class was calculated for the Resistomap versus CARD/RGI at percentages of 25%, 50%, and 95%.

Taxonomic classification
Taxonomic classification of the reads of the four samples was labelled using the kraken tool (Fig 1), the distribution of taxa shows that the most abundant bacterial groups in all four samples are Gammaproteobacteria, Actinomycetes, Betaproteobacteria, and Alphaproteobacteria, in that order.However, their relative abundances differ between samples (S1 Table ).For example, Sample 1 (S1) and Sample 2 (S2) have similar proportions of Gammaproteobacteria, while Sample 3 (S3) and Sample 4 (S4) have lower proportions.Compared to the other three samples, S4 shows a relatively lower proportion of Gammaproteobacteria, Betaproteobacteria, and Epsilonproteobacteria, while it has a higher proportion of Alphaproteobacteria, Actinomycetes, and Flavobacteriia.S4 also has a higher proportion of Acidimicrobiia, Thermoleophilia, Deinococci, and Planctomycetia compared to the other three samples.Overall, the microbial community in S4 appears to differ significantly from the other three samples and seems to constitute a more diverse microbial community.A table that provides data on the abundance of each of these bacterial groups for each sample is provided as supplementary data (S2 Table ).
In order to ensure the comprehensive coverage of low abundance species/genes, we integrated rarefaction into the CARD/RGI workflow (Fig 2).Leveling off of rarefaction curves was observed, in particular for drug classes (Fig 2 ), indicating that the shotgun sequencing effort is likely nearly comprehensive.However, the rarefaction curves also showed strong sensitivity to percentage-coverage thresholding.For example, far more AMR genes are detected for 10% percentage coverage than for 100% percentage coverage.For this reason, we determined that the percentage coverage parameter was critical in determining if an ARG was present or absent in a sample.

Gene ambiguity as reported by Resistomap and curated in the CARD database
In order to compare the shotgun sequencing method with the qPCR method, we needed to determine what genes in the CARD databases matched with the primer pairs used by Resistomap.This is also necessary as genes tend to have multiple names (synonyms) when described by different research groups.When we blasted the Resistomap primers pairs against the CARD database, we found that 93 primer pairs matched a unique gene, while 45 primer pairs matched multiple genes in the database.In the latter case, these multiple targets are of the same or similar function, mostly paralogs (S2 Table, see all_multiTargets.csv in GitHub repository).In 37 of these 45 multiple-match primer pairs, all the identified genes had a common ARO term as a parent.For the other seven primers, a common grandparent was identified.For example, Resistomap's primer '16S new 2' targets both Mycobacterium tuberculosis 16S rRNA mutation conferring resistance to capreomycin (ARO:3003411) and Escherichia coli 16S rRNA (rrnB) mutation conferring resistance to tetracycline' (ARO:3003211).In one case, two matched genes (CMY-8b and MOX14) had only a common great-grandparent in ARO, but both confer resistance to beta-lactam antimicrobials.We limited our comparison to those cases, where primers matched one or more CARD gene targets, henceforth referred to as single-target and multi-target, respectively.Tables that detail all the genes that were identified by CARD/RGI but not Resistomap, as well as genes identified by Resistomap but not CARD/RGI, are provided in the supplementary data (S3 and S4 Tables). .Supplementary data that indicates the top genes identified in each class using CARD/RGI analysis is provided (S5 Table ).ARGs that confer resistance to all 11 antibiotic classes were detected in the 4 samples using both methods.The Resistomap qPCR method identified an average of 115 genes for samples S1, S2, and S3, with a visibly lower number of ARGs (94) observed in sample S4.The low number of ARGs detected in sample S4 was particularly pronounced when the CARD/RGI method was used as only 120 ARGs were detected in sample S4 compared to an average of 342 genes for S1, S2, and S3.The most abundant ARG classes conferred resistance to aminoglycoside, MDR, MLSB, tetracycline and beta-lactams.The difference in the abundance of detected ARGs across the different tested locations was more apparent using the CARD/RGI method.For example, a similar number of ARGs against Beta-lactam were detected in all 4 samples using Resistomap's qPCR-based method, while the number of beta-lactam ARGs detected with the CARD/RGI method varied significantly (S1 = 53, S2 = 35, S3 = 95 and S4 = 9).Notably, when considering the genes that were identified by CARD/RGI but not by Resistomap 46 of those genes belonged to the beta-lactamase family, some of which (e,g.NPS-1 and LCR-1), were identified with high abundance values of up to 247.89 and 132.76, respectively, in S1, as well as OXA-732, which reached up to 201.87 in S3.Vancomycin resistance ARGs represented the least abundant class of ARGs.Indeed, the Resistomap qPCR method identified only one ARG in sample S1, two ARGs in sample S2, and none in samples S3 and S4.Very similarly, the CARD/RGI method identified one ARG in sample S1, one ARG in sample S2, two ARGs in sample S3, and none observed in sample S4.Considering how important vancomycin is as one of the last resort antibiotics, it would be very important to monitor trends in vancomycin ARGs in wastewater, even if they are not very abundant.
The abundance of ARGs that were identified by both methods was used to plot a heatmap (Fig 4) in order to further compare these methods.The Resistomap PCR primer pairs were matched to the CARD genes (see Methods and Section 3.3 above) in order to allow for this comparison.The heatmap shows and compares the relative abundances of 93 genes in each of the tested samples obtained using Resistomap qPCR or sequencing and CARD/RGI, represented as normalized gene coverage, where detection coverage thresholds for each ARG are set at either 95% or 50%.
A method-dependent difference in the number of detected genes was observed for all tested samples.In the case of S1, 23 out of 93 genes were detected by Resistomap while 24 and 33 genes were detected by CARD/RGI >95% and CARD/RGI >50%, respectively.A higher number of genes was also detected by Resistomap (28 genes) and CARD/RGI >50% (30 genes) compared to CARD/RGI >95% (19 genes) in S2.In contrast, more genes were detected by Resistomap in S3 and S4 samples.The CARD/RGI results for gene read percentage coverages of 95% and 50% showed a very similar profile with differences in the cases of genes showing low abundance and only being detected at a percentage coverage of 50%.One exception was mexE gene (MDR resistance) which was detected only with a gene read percentage coverage of 50% and whose relative abundance ranged between 2.5 and 41 in the tested samples.Nine genes displayed similar abundance patterns and fold coverage values using both methods across all samples (tetM, mefA, strB, aadB in CARD/RGI 95% and tetM, mefA, InuC, strB, aadB, aac(6')_II, acrF, oqxA, and mdtA in CARD/RGI 50%), and 49 genes were not detected by any of the methods including all ARGs conferring resistance to vancomycin (Table 1).Among the 93 genes, 8 genes were exclusively detected by Resistomap (aadD, apmA, catP, adeI, dfrA27, dfrA10, tetE and blaTLA) and showed low relative abundance, while 6 genes were only identified by CARD/RGI (tetS, ermA/ermTR, aacA/aphD, addA7, catB3 and dfrG) with a fold coverage between 3-74.Genes showing high abundance (�0.02 and �100 in Resistomap and CARD/RGI methods, respectively) included tetM (in S2) which confers resistance to tetracycline, mefA (in S1 and S2) which confers resistance to MLSB, aadB (in S3) which confers resistance to aminoglycoside, and catQ (in S2) which confers resistance to phenicol.
To understand why the mexB gene was only identified by sequencing but not qPCR, a mapping of the gene reads against the CARD gene reference was performed (Fig 5).
The mapping results showed uneven fold coverage of the mexB gene.By further investigating the gene sequence targeted by Resistomap primers, a large number of mutations (nongrey) is observed with three mutations occurring at the first three positions targeted by the reverse primer (Fig 6).In addition, we identified 45 primer pairs that matched with multiple CARD targets (multitarget), however, these targets are homologous (or are closely related in terms of the ARO graph distance, see methods) (Fig 7).The number of detected genes in each sample was comparable between Resistomap and CARD/RGI (50%) methods while CARD/RGI (95%) showed lower numbers.Four genes were exclusively detected using the CARD/RGI method with a fold coverage between 0.52-87, and 7 genes were only detected by Resistomap (<0.01 relative abundance) (Table 1).Genes tetW (tetracycline), aadA16 and aadA6 (aminoglycoside) were abundant in samples S1, S2 and S3.Gene blaGES (beta-lactam) was also abundant in sample S3.
The results showed a reasonable correlation in the identification of genes by both methods.The unique primer matches scored a Pearson correlation coefficient of 0.63 while the multitarget primers scored 0.8.The breakdown of the correlation by antibiotic class (Fig 8 ) showed a strong correlation (> 0.7) between the results of the two methods for tetracycline, MLSB, aminoglycoside, phenicol, quinolone, and beta-lactam.In addition, both methods agree on vancomycin not being detected.In contrast, the MDR class showed a negative correlation while trimethoprim and others showed negative correlations only in the Resistomap vs CARD 95% comparison.

Discussion
Quantitative techniques such as qPCR and sequencing (next-generation and nanopore) are being increasingly used to monitor the presence of ARGs in environmental samples and identify new resistance markers [39][40][41][42].While qPCR represent a fast, cost-effective and sensitive approach to detect and quantify ARGs, the use of this technique is dependent upon known target sequences that contain conserved primer target sites.However, the mentioned preconditions are not always satisfied and shotgun metagenomic sequencing is sometimes required for  Identified ARGs by both methods conferred resistance mainly to aminoglycoside, MDR, MLSB, tetracycline and beta-lactams, which are often detected in wastewater and are used extensively in human or veterinary medicine, including in-feed [43,44].Among the detected ARGs, some pose a serious threat to public health.The detection of mefA ARG in all of the tested samples at relatively high abundance is clinically alarming.The gene confers resistance  [45,46].Antibiotics are more frequently used in healthcare facilities making the ARGs profile in the wastewater of hospitals richer compared to other locations [47].This was reflected by more genes showing a higher abundance in sample S3 (healthcare facility) compared to the other samples regardless of the used analysis method.
More variation in the distribution of genes among samples was observed by non-targeted sequencing highlighting the differences in the local environments of the tested samples and demonstrating that sequencing can provide more information on the ARGs composition in wastewater compared to HT-qPCR.ARGs conferring resistance to beta-lactams for example were more prevalent in the wastewater sample collected nearby a healthcare facility by sequencing this sample.Beta-lactam antibiotics have a broad spectrum of activity and are frequently prescribed in clinical settings [48].Their extensive use raises the likelihood of the emergence of resistant strains as observed by the increasing trend of percentage resistant isolates in the UAE during the period 2010-2020 [49].It has been previously found that up to 37% of hospital wastewater samples contained ARGs conferring resistance to beta-lactams compared to 18% in municipal wastewater [50].
Another ARG class showing significant method-dependent variation is vancomycin.Overall, genes conferring resistance to this antibiotic were the least detected across all samples.Vancomycin inhibits cell-wall synthesis in Gram-positive bacteria and is commonly used against streptococci and enterococci resistant to beta-lactams and penicillin [51].The development of vancomycin-resistant strains is rare with very few known examples including vancomycinresistant enterococci (Enterococcus faecalis and Enterococcus faecium) and vancomycin-resistant Staphylococcus aureus, explaining the low number of detected ARGs [52,53].This was further supported by the report showing that Enterococcus spp.and S. aureus isolates from different sources in the UAE showed 98 and 100% susceptibility to vancomycin, respectively [49].These organisms are considered a leading cause of nosocomial infections, which possibly explains the detection of vancomycin ARGs in sample S3 by sequencing [54].However, vancomycin ARGs were also detected in samples S1 and S2 by Resistomap HT-qPCR and further analysis is required to evaluate the significance of this result from a public health perspective.
A possible explanation for the observed discrepancies in the results obtained by the two different methods is that by sequencing, the reads obtained are aligned to all ARGs-like sequences in CARD database.In contrast, HT-qPCR detection is limited by the used primers which are designed to target sections of ARG sequences.In addition, HT-qPCR is more vulnerable to the presence of inhibitors in the samples and therefore the used condition might not be optimal to all ARGs primers resulting in lower abundances than the actual values in wastewater.By analyzing a large number of wastewater [32], and fecal samples [55,56] by qPCR, potential issues regarding the specificity of used primer sets have been raised.
A strong positive correlation was obtained between the relative abundance (obtained by Resistomap HT-qPCR) and the fold-coverage (obtained by TruSeq DNA sequencing with 50% percentage coverage) of all detected ARG classes except for MDR.Three ARGs are responsible for the observed negative correlation for this class of antibiotics.The ARG oprD was detected by Resitomap in samples S1-S3 while only in S3 by sequencing.HT-qPCR also exclusively detected the gene adeI at low abundance in sample S3.On the other hand, the abundance of gene mexE in samples S1-S3 was significantly higher by sequencing and gene mexB was almost exclusively detected at high abundance by CARD/RGI.P. aeruginosa with mutations in mexE and mexB overexpress transport-proteins involved in the excretion of molecules from the bacterial cell and result in reducing antibiotic susceptibility [57][58][59].Antibiotic-resistant P. aeruginosa is another leading causative agent for nosocomial infection including wound and urinary infections, and ventilator-associated pneumonia [60], and it was found at high abundance in sample S3, collected nearby a healthcare facility.However, this organism is highly adaptive and can colonize a broad range of environments [61].Mapping the obtained gene reads of mexB against the CARD gene reference revealed uneven fold coverage that can be explained by the regions with higher coverage being common among different species.Fold-coverage is therefore an overestimate and an ideal algorithm should differentiate fold-coverage for specific and unspecific regions.In addition, the used primers can be ineffective by selectively targeting taxa-specific ARG and therefore under detecting the resistance marker in the sample.The mexB gene also contained a large number of mutations that could have affected the specificity of the used Resistomap primer, especially the three mutations in the region targeted by the reverse primer (Fig 6).
It was notable that Staphylococci were not detected in any of the four water samples using either method (Table 1).In one previous study that examined 62 wastewater samples, Staphylococcus was recovered from only 12 (19%) of the samples [62].Another study utilizing qPCR based methods similarly found that S. aureus is not very prevalent in wastewater [63].In our own study where 88 wastewater samples were analyzed, Staphylococcus was only detected in 2 out of the 88 samples (2.3%) [29].Furthermore, the abundance of Staphylococcus in those two samples was very low (0.03% and 0.06%) of the bacterial community.Taking these data together, the lack of detection could indicate that low abundance microbes in wastewater require either more sequencing depth or enrichment prior to DNA extraction if they are an important component of the AMR surveillance program.
Upon comparison of the Resistomap and CARD/RGI methods, it can be concluded that qPCR has the advantage of having perfect matches in a large number of CARD genes (2787) despite using only 384 targets.In addition, PCR-based methods like Resistomap are highly sensitive to low abundance genes, but are vulnerable to mutations in primer target region (introducing false negatives).Another disadvantage of Resistomap is that it identifies genes by only matching the primers' sequences, often only covering a short target region, thus not detecting partial genes and potential loss of function.On the other hand, CARD/RGI contains 4891 genes but is more likely to miss genes due to only partial coverage resulting from large microbial complexity.For this reason, it is important to set a lower threshold at first and analyze genes of interest to determine if they are functional or not depending on coverage, a function that is possible using bioinformatics pipelines.Therefore, each method has its advantages and disadvantages, and the choice of the appropriate method depends on the specific research objectives and available resources.
The challenges presented by this study include the issue of conflicting matches between different platforms when using different antimicrobial gene naming schemes.This highlights the need for standardization in naming conventions across different databases and platforms.Additionally, the limitation of CARD/RGI in ignoring important antimicrobial resistance mechanisms such as MGE is a significant challenge for researchers and clinicians.This highlights the importance of continued development and improvement of these databases to encompass a broader range of resistance mechanisms.The complexity of the relationship between primers and genes in databases like CARD is another challenge that requires careful consideration to ensure accurate identification and interpretation of resistance genes.Overall, these challenges highlight the need for ongoing research and development to enhance the accuracy and utility of antimicrobial resistance databases for clinical and research applications.In addition, these challenges also indicate the utilization of both methods may be deemed inappropriate in environments lacking essential bioinformatics proficiency.Even if samples analysis are outsourced, it would be advisable to hire individuals with expertise within municipal services or entities that can understand these reports.This last point reinforces the need to train future microbiologists in bioinformatics/statistics and for municipalities to ensure they budget for the infrastructure and workforce needed to monitor wastewater for AMR successfully.
In conclusion, our study highlights the importance of using both PCR and NGS methods for gene detection and identification.While both methods identify some genes with similar abundance, they also complement each other by identifying genes that are missed by the other method.The combination of both methods can provide more accurate and comprehensive results, and careful analysis and optimization of the threshold settings can help improve the reliability of the results.Overall, considering the strengths and limitations of each method and using those in a complementary manner is essential for optimal results.

Fig 1 .
Fig 1.The distribution of reads assigned to "Class" taxon using the Kraken tool for the four wastewater samples.S1 was collected from an urban residential site, S2 from an industrial site, S3 from a healthcare facility, and S4 was collected from a rural sewage treatment plant influent.https://doi.org/10.1371/journal.pone.0298325.g001

Fig 2 .Fig 3 .
Fig 2. Rarefaction landscape using different gene coverage thresholds for number of antibiotic resistance ontology (ARO) terms, drug classes and AMR gene families for sample S1.Each line graph contains lines representing the different gene cover thresholds used in the rarefaction analysis (indicated in the legend).https://doi.org/10.1371/journal.pone.0298325.g002

Fig 5 .Fig 6 .
Fig 5. CARD/RGI alignment results showing reads coverage of gene mexB.Note that more than 95% of the gene is covered with at least one read in all samples.https://doi.org/10.1371/journal.pone.0298325.g005

Table 1 .
(Continued) # mexB was only detected by CARD/RGI in certain samples with a high relative abundance value.It was therefore included under only CARD/RGI.Only CARD/RGI column includes genes with percentage-coverages of 95% and 50%.https://doi.org/10.1371/journal.pone.0298325.t001